Example 4: Currency Substitution with capital accumulation. Two sector economy with productivity shock in traded sector. See paper for details.
Manoj Atolia: September 18, 2007.
Works with both Mathematica 6.0 and 5.2.
Trial Version

Off[General::spell];Off[General::spell1];Off[Unset::norep];Off[FindMinimum::lstol];Off[FindMinimum::reged];Off[InterpolatingFunction::dmval];Off[Precision::precsm];

(* The precision of preliminary calculations. *)
{precision = 100, $MinPrecision=digits=precision};

Initialization

(* Set the directory where Mathematica will store intermediate results. *)
SetDirectory["c:\Documents and Settings\Manoj Atolia\My Documents\Academic\Research\CurrentResearch\Shooting Methods\Computations\Reverse Shooting\Recent Version Files of Examples - Package Form"];

(* ALWAYS ENTER NUMERICAL VALUES AS FRACTIONS AND NOT AS DECIMALS TO ENSURE REQUISITE PRECISION IN CALCULATIONS *)   
{τ=(1/2),ρ=(1/10),σ=2,δ=(1/20),β=(1/2),γ=(1/2),α=(1/2),Ω=2}//N                                     (* Parameters of the model *)

{(Subscript[bbyY, o])=(25/100),(Subscript[mbyF, o])=(1/2),(Subscript[mbyY, o])=(15/200),(Subscript[π, o])=(1/2),(Subscript[θn, o])=(1/2),(Subscript[θt, o])=(1/2),Pt=1}//N           (* Data for calibration *)

{(Subscript[z, o])=1,lbar=1,(Subscript[Pn, o])=1}//N                                                                                                    (* Normalization *)

{(Subscript[z, f])=(6/5)(Subscript[z, o]),ω=(1/5)}//N                                                                                            (* Shock - Exogenous rise in productivity *)

Calibration and Steady State

(* Define utility function, production function etc. *)
V[Ex_,Pn_]:=(1/(1-(1/τ)))((((Ex/((γ^γ)(((1-γ))^(1-γ))(Pn^γ)))))^(1-(1/τ)));(Subscript[V, 1])[Ex_,Pn_]:=Evaluate[D[V[Ex,Pn],Ex]];(Subscript[V, 2])[Ex_,Pn_]:=Evaluate[D[V[Ex,Pn],Pn]];
ϕ[vm_,vF_]:=((Subscript[g, 0])/(1-1/τ))(((((((Subscript[g, 1])(((vm))^((σ-1)/σ))+(1-(Subscript[g, 1]))(((vF))^((σ-1)/σ))))^(σ/(σ-1)))))^(1-1/τ));(Subscript[ϕ, 1])[vm_,vF_]:=Evaluate[D[ϕ[vm,vF],vm]];(Subscript[ϕ, 2])[vm_,vF_]=Evaluate[D[ϕ[vm,vF],vF]];
F[k_,l_]:=((((Subscript[an, 0])(l^((β-1)/β))+(Subscript[an, 1])(k^((β-1)/β))))^(β/(β-1)));(Subscript[F, 1])[k_,l_]:=Evaluate[D[F[k,l],k]];(Subscript[F, 2])[k_,l_]:=Evaluate[D[F[k,l],l]];
G[k_,l_]:=((((Subscript[at, 0])(l^((β-1)/β))+(Subscript[at, 1])(k^((β-1)/β))))^(β/(β-1)));(Subscript[G, 1])[k_,l_]:=Evaluate[D[G[k,l],k]];(Subscript[G, 2])[k_,l_]:=Evaluate[D[G[k,l],l]];
ψ[i_,k_]:=((Subscript[a, 2])/2)((((i/k)-δ))^2);(Subscript[a, 2])=(1/(Ω δ));ψ'[i_,k_]:=(Subscript[a, 2])(i/k-δ);
ψ''[i_,k_]:=(Subscript[a, 2]);
Pk[Pn_]:=(α^α)(((1-α))^(1-α))(Pt^(1-α))(Pn^α);

(* Calibrate and compute the original steady state *)
{(Subscript[Yn, o])=(Subscript[z, o])F[(Subscript[kn, o]),(Subscript[ln, o])];(Subscript[Yt, o])=(Subscript[z, o])G[(Subscript[kt, o]),lbar-(Subscript[ln, o])];(Subscript[Y, o])=(Subscript[Pn, o])(Subscript[Yn, o])+Pt (Subscript[Yt, o]);(Subscript[Ex, o])=(Subscript[Y, o])-Pk[(Subscript[Pn, o])]δ ((Subscript[kn, o])+(Subscript[kt, o])),(Subscript[m, o])= (Subscript[mbyY, o]) (Subscript[Y, o]),(Subscript[F, o])=((Subscript[m, o])/(Subscript[mbyF, o])),(Subscript[b, o])=(Subscript[bbyY, o]) (Subscript[Y, o]),L=(Subscript[π, o])(Subscript[m, o])-ρ (Subscript[b, o]),(Subscript[v, o])=((Subscript[Pn, o])^(-γ))};(Subscript[an, 0])=1;
{(Subscript[an, 1]),(Subscript[at, 0]),(Subscript[at, 1]),(Subscript[ln, o]),(Subscript[kn, o]),(Subscript[kt, o]),(Subscript[Ex, o])}=.;({(Subscript[an, 1]),(Subscript[at, 0]),(Subscript[at, 1]),(Subscript[ln, o]),(Subscript[kn, o]),(Subscript[kt, o]),(Subscript[Ex, o])}={(Subscript[an, 1]),(Subscript[at, 0]),(Subscript[at, 1]),(Subscript[ln, o]),(Subscript[kn, o]),(Subscript[kt, o]),(Subscript[Ex, o])}/.FindRoot[{(Subscript[z, o])(Subscript[Pn, o])(Subscript[F, 1])[(Subscript[kn, o]),(Subscript[ln, o])] == (ρ+δ)Pk[(Subscript[Pn, o])],(Subscript[θn, o])==(((Subscript[z, o])(Subscript[F, 1])[(Subscript[kn, o]),(Subscript[ln, o])](Subscript[kn, o]))/((Subscript[z, o])F[(Subscript[kn, o]),(Subscript[ln, o])])),(Subscript[z, o])(Subscript[G, 1])[(Subscript[kt, o]),lbar-(Subscript[ln, o])] == (ρ+δ)Pk[(Subscript[Pn, o])],(Subscript[θt, o])==(((Subscript[z, o])(Subscript[G, 1])[(Subscript[kt, o]),lbar-(Subscript[ln, o])](Subscript[kt, o]))/((Subscript[z, o])G[(Subscript[kt, o]),lbar-(Subscript[ln, o])])),(Subscript[Pn, o])(Subscript[z, o])(Subscript[F, 2])[(Subscript[kn, o]),(Subscript[ln, o])]==(Subscript[z, o])(Subscript[G, 2])[(Subscript[kt, o]),lbar-(Subscript[ln, o])],(Subscript[z, o])F[(Subscript[kn, o]),(Subscript[ln, o])]==((γ (Subscript[Ex, o]))/(Subscript[Pn, o]))+((α Pk[(Subscript[Pn, o])])/(Subscript[Pn, o]))δ((Subscript[kn, o])+(Subscript[kt, o])),(Subscript[Ex, o])==(Subscript[Y, o])-Pk[(Subscript[Pn, o])]δ ((Subscript[kn, o])+(Subscript[kt, o]))},{(Subscript[an, 1]),1},{(Subscript[at, 0]),1},{(Subscript[at, 1]),1},{(Subscript[ln, o]),1/2},{(Subscript[kn, o]),1},{(Subscript[kt, o]),1},{(Subscript[Ex, o]),1/2},MaxIterations->1000,AccuracyGoal->digits,WorkingPrecision->digits])//N
{(Subscript[g, 0])=.,(Subscript[g, 1])=.};eq1:= (Subscript[v, o])(Subscript[ϕ, 1])[(Subscript[v, o])(Subscript[m, o]),(Subscript[v, o])(Subscript[F, o])]== (Subscript[V, 1])[(Subscript[Ex, o]),(Subscript[Pn, o])](ρ+(Subscript[π, o]));eq2:=(Subscript[v, o])(Subscript[ϕ, 2])[(Subscript[v, o])(Subscript[m, o]),(Subscript[v, o])(Subscript[F, o])]== (Subscript[V, 1])[(Subscript[Ex, o]),(Subscript[Pn, o])]ρ ;({(Subscript[g, 0]),(Subscript[g, 1])}={(Subscript[g, 0]),(Subscript[g, 1])} /. FindRoot[{eq1,eq2},{(Subscript[g, 0]),(10^(-6))},{(Subscript[g, 1]),1},MaxIterations->1000,AccuracyGoal->$MinPrecision,WorkingPrecision->digits])//N

(* Compute the new steady state *)
{(Subscript[Yn, f])=(Subscript[z, o])F[(Subscript[kn, f]),(Subscript[ln, f])];(Subscript[Yt, f])=(Subscript[z, f])G[(Subscript[kt, f]),lbar-(Subscript[ln, f])];(Subscript[Y, f])=(Subscript[Pn, f])(Subscript[Yn, f])+Pt (Subscript[Yt, f]);(Subscript[Ex, f])=(Subscript[Y, f])-Pk[(Subscript[Pn, f])]δ ((Subscript[kn, f])+(Subscript[kt, f]));(Subscript[b, o])=(Subscript[bbyY, o]) (Subscript[Y, o])};(Subscript[v, f])=((Subscript[Pn, f])^(-γ));
{(Subscript[Pn, f]),(Subscript[ln, f]),(Subscript[kn, f]),(Subscript[kt, f]),(Subscript[Ex, f])}=.;({(Subscript[Pn, f]),(Subscript[ln, f]),(Subscript[kn, f]),(Subscript[kt, f]),(Subscript[Ex, f])}={(Subscript[Pn, f]),(Subscript[ln, f]),(Subscript[kn, f]),(Subscript[kt, f]),(Subscript[Ex, f])}/.FindRoot[{(Subscript[z, o])(Subscript[Pn, f])(Subscript[F, 1])[(Subscript[kn, f]),(Subscript[ln, f])] == (ρ+δ)Pk[(Subscript[Pn, f])],(Subscript[z, f])(Subscript[G, 1])[(Subscript[kt, f]),lbar-(Subscript[ln, f])] == (ρ+δ)Pk[(Subscript[Pn, f])],(Subscript[Pn, f])(Subscript[z, o])(Subscript[F, 2])[(Subscript[kn, f]),(Subscript[ln, f])]==(Subscript[z, f])(Subscript[G, 2])[(Subscript[kt, f]),lbar-(Subscript[ln, f])],(Subscript[z, o])F[(Subscript[kn, f]),(Subscript[ln, f])]==((γ (Subscript[Ex, f]))/(Subscript[Pn, f]))+((α Pk[(Subscript[Pn, f])])/(Subscript[Pn, f]))δ((Subscript[kn, f])+(Subscript[kt, f])),(Subscript[Ex, f])==(Subscript[Y, f])-Pk[(Subscript[Pn, f])]δ ((Subscript[kn, f])+(Subscript[kt, f]))},{(Subscript[Pn, f]),1},{(Subscript[ln, f]),1/2},{(Subscript[kn, f]),1},{(Subscript[kt, f]),1},{(Subscript[Ex, f]),1/2},MaxIterations->1000,AccuracyGoal->digits,WorkingPrecision->digits])//N
{(Subscript[m, f]),(Subscript[F, f]),(Subscript[π, f])}=.;{(Subscript[m, f]),(Subscript[F, f]),(Subscript[π, f])}={(Subscript[m, f]),(Subscript[F, f]),(Subscript[π, f])}/.FindRoot[{(((Subscript[v, f])(Subscript[ϕ, 2])[(Subscript[v, f])(Subscript[m, f]),(Subscript[v, f])(Subscript[F, f])])/((Subscript[V, 1])[(Subscript[Ex, f]),(Subscript[Pn, f])]))==ρ,(((Subscript[v, f])(Subscript[ϕ, 1])[(Subscript[v, f])(Subscript[m, f]),(Subscript[v, f])(Subscript[F, f])])/((Subscript[V, 1])[(Subscript[Ex, f]),(Subscript[Pn, f])]))==ρ+(Subscript[π, f]),(Subscript[π, f])(Subscript[m, f])==L+ρ (Subscript[b, o]) },{(Subscript[m, f]),(Subscript[m, o])},{(Subscript[π, f]),(Subscript[π, o])},{(Subscript[F, f]),(Subscript[F, o])},AccuracyGoal->$MinPrecision,WorkingPrecision->digits];
{{(Subscript[m, o]),(Subscript[F, o]),(Subscript[kn, o]),(Subscript[kt, o]),(Subscript[Ex, o]),(Subscript[Pn, o]),(Subscript[ln, o]),(Subscript[Y, o]),(Subscript[b, o]),L,(Subscript[z, o]),(Subscript[π, o]),(Subscript[Ex, o])/(Subscript[m, o]),(Subscript[m, o])/(Subscript[Y, o]),(Subscript[m, o])/(Subscript[F, o])},{(Subscript[m, f]),(Subscript[F, f]),(Subscript[kn, f]),(Subscript[kt, f]),(Subscript[Ex, f]),(Subscript[Pn, f]),(Subscript[ln, f]),(Subscript[Y, f]),(Subscript[b, o]),L,(Subscript[z, f]),(Subscript[π, f]),(Subscript[Ex, f])/(Subscript[m, f]),(Subscript[m, f])/(Subscript[Y, f]),(Subscript[m, f])/(Subscript[F, f])}}//N//MatrixForm

Global Nonlinear Saddlepath

(* Solve for the endogenous variables *)
in[Pn_,λ1_,λ2_,kn_]=Evaluate[in/.Solve[λ2==λ1 Pk[Pn](1+ψ'[in,kn]),in][[1]]];
it[Pn_,λ1_,λ3_,kt_]=Evaluate[it/.Solve[λ3==λ1 Pk[Pn](1+ψ'[it,kt]),it][[1]]];
i[Pn_,λ1_,m_,F_]=(((Pn^(-γ))(Subscript[ϕ, 1])[(Pn^(-γ))m,(Pn^(-γ))F])/(λ1));
χ[Pn_,λ1_,m_,F_]=(((Pn^(-γ))((Subscript[ϕ, 1])[(Pn^(-γ))m,(Pn^(-γ))F]-(Subscript[ϕ, 2])[(Pn^(-γ))m,(Pn^(-γ))F]))/(λ1));
Ex[Pn_,λ1_]=Evaluate[Ex/.Solve[(Subscript[V, 1])[Ex,Pn]==λ1,Ex]][[2]];
{{(Subscript[λ1, f])=(Subscript[V, 1])[(Subscript[Ex, f]),(Subscript[Pn, f])],(Subscript[λ1, o])=(Subscript[V, 1])[(Subscript[Ex, o]),(Subscript[Pn, o])]},{{(Subscript[λ2, f]),(Subscript[λ3, f])}={(Subscript[λ1, f])Pk[(Subscript[Pn, f])],(Subscript[λ1, f])Pk[(Subscript[Pn, f])]},{(Subscript[λ2, o]),(Subscript[λ3, o])}={(Subscript[λ1, o])Pk[(Subscript[Pn, o])],(Subscript[λ1, o])Pk[(Subscript[Pn, o])]}}}//N

(* Specify the nonlinear system with state variables at the end and maintain same order subsequently *)
λ1dot = (ρ-(i-χ))λ1;
λ2dot = (ρ+δ)λ2 - λ1(Pn (Subscript[F, 1])[kn,ln]-Pk[Pn](ψ[in,kn]-ψ'[in,kn]in/kn));
λ3dot=(ρ+δ)λ3-λ1(z  (Subscript[G, 1])[kt,lbar-ln]-Pk[Pn](ψ[it,kt]-ψ'[it,kt]it/kt));
mdot=L+(i-χ) (Subscript[b, o])-χ m;
Fdot = Pn F[kn, ln] + z G[kt, lbar - ln]-Ex - Pk[Pn](in + ψ[in,kn]kn + it + ψ[it, kt]kt);
kndot=in-δ kn;
ktdot=it-δ kt;
zdot =  ω((Subscript[z, f])-z);
eq1:=Pn (Subscript[F, 2])[kn,ln]==z (Subscript[G, 2])[kt,lbar-ln];
eq2:=F[kn,ln]==((γ Ex[Pn,λ1])/Pn)+((α Pk[Pn])/Pn)(in[Pn,λ1,λ2,kn] +ψ[in[Pn,λ1,λ2,kn],kn]kn+it[Pn,λ1,λ3,kt]+ψ[it[Pn,λ1,λ3,kt],kt]kt);
{Pndot,lndot}={Dt[Pn],Dt[ln]}/.Solve[{Dt[eq1]/.{Dt[λ1]->λ1dot,Dt[λ2]->λ2dot,Dt[λ3]->λ3dot,Dt[m]->mdot,Dt[F]->Fdot,Dt[kn]->kndot,Dt[kt]->ktdot,Dt[z]->zdot},Dt[eq2]/.{Dt[λ1]->λ1dot,Dt[λ2]->λ2dot,Dt[λ3]->λ3dot,Dt[m]->mdot,Dt[F]->Fdot,Dt[kn]->kndot,Dt[kt]->ktdot,Dt[z]->zdot}},{Dt[Pn],Dt[ln]}][[1]];
system={Pndot, lndot, λ1dot, λ2dot, λ3dot, mdot, Fdot,kndot,ktdot,zdot};

(* INITIALIZATION CELL - Initialization for finding the required trajectory. *)
{   T = 20;          (* solution time (years). BEWARE!!! Setting Tss larger than necessary can substantially increase computational time.*) 
precision = 16;        (* starting value of precision of calculations*) ;
delta =  (10^(-1))(((((((Subscript[F, f])-(Subscript[F, o])))^2)+((((Subscript[kn, f])-(Subscript[kn, o])))^2)+((((Subscript[kt, f])-(Subscript[kt, o])))^2)+((((Subscript[z, f])-(Subscript[z, o])))^2)))^(1/2));(* size of δ-simplex  or distance of simplex  from the final steady state *);
tol = (10^(-3));                         (* desired accurcy of the solutions of the variables in the system *);};

?solve

(* MAIN CELL - Solving for the global nonlinear saddlepath. *)
Timing[{gtmin,gminsol}=solve[system,{t},{Pn,ln,λ1,λ2,λ3,m,F,kn,kt,z},{(Subscript[Pn, o]),(Subscript[ln, o]),(Subscript[λ1, o]),(Subscript[λ2, o]),(Subscript[λ3, o]),(Subscript[m, o]),(Subscript[F, o]),(Subscript[kn, o]),(Subscript[kt, o]),(Subscript[z, o])},{(Subscript[Pn, f]),(Subscript[ln, f]),(Subscript[λ1, f]),(Subscript[λ2, f]),(Subscript[λ3, f]),(Subscript[m, f]),(Subscript[F, f]),(Subscript[kn, f]),(Subscript[kt, f]),(Subscript[z, f])},{{in,{Pn,λ1,λ2,kn}},{it,{Pn, λ1, λ3, kt}},{i,{Pn, λ1, m, F}},{χ,{Pn,λ1, m, F}},{Ex,{Pn,λ1}}},{T,tol,precision,delta}]]

(* Plot the final solutions. *)
tplot = 31;
mPlot=Plot[{If[t<0,(Subscript[m, o]),m[gtmin-t]]/.gminsol,(Subscript[m, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",m},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
ExPlot=Plot[{If[t<0,(Subscript[Ex, o]),Ex[Pn[gtmin-t],λ1[gtmin-t]]]/.gminsol,(Subscript[Ex, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",Ex},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
FPlot=Plot[{If[t<0,(Subscript[F, o]),F[gtmin-t]]/.gminsol,(Subscript[F, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",F},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
inPlot=Plot[{If[t<0,δ (Subscript[kn, o]),in[Pn[gtmin-t],λ1[gtmin-t],λ2[gtmin-t],kn[gtmin-t]]]/.gminsol,δ (Subscript[kn, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",(Subscript[I, n])},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
knPlot=Plot[{If[t<0,(Subscript[kn, o]),kn[gtmin-t]]/.gminsol,(Subscript[kn, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",(Subscript[k, n])},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
lnPlot=Plot[{If[t<0,(Subscript[ln, o]),ln[gtmin-t]]/.gminsol,(Subscript[ln, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",(Subscript[l, n])},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
itPlot=Plot[{If[t<0,δ (Subscript[kt, o]),it[Pn[gtmin-t],λ1[gtmin-t],λ3[gtmin-t],kt[gtmin-t]]]/.gminsol,δ (Subscript[kt, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",(Subscript[I, t])},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
ktPlot=Plot[{If[t<0,(Subscript[kt, o]),kt[gtmin-t]]/.gminsol,(Subscript[kt, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",(Subscript[k, t])},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
YPlot=Plot[{If[t<0,(Subscript[Y, o]),Pn[gtmin-t]F[kn[gtmin-t],ln[gtmin-t]]+z[gtmin-t]G[kt[gtmin-t],lbar-ln[gtmin-t]]]/.gminsol,(Subscript[Y, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",Y},AxesOrigin->{-tplot/10,Automatic},
DisplayFunction->Identity];
iPlot=Plot[{If[t<0,ρ+(Subscript[π, o]),i[Pn[gtmin-t],λ1[gtmin-t],m[gtmin-t],F[gtmin-t]]]/.gminsol,ρ+(Subscript[π, f])},{t,-tplot/10,tplot},PlotRange->All,AxesOrigin->{-tplot/10,Automatic},AxesLabel->{"",i},DisplayFunction->Identity];
χPlot=Plot[{If[t<0,(Subscript[π, o]),χ[Pn[gtmin-t],λ1[gtmin-t],m[gtmin-t],F[gtmin-t]]]/.gminsol,(Subscript[π, f])},{t,-tplot/10,tplot},PlotRange->All,AxesOrigin->{-tplot/10,(Subscript[π, f])},AxesLabel->{"",χ},DisplayFunction->Identity];zPlot=Plot[{If[t<0,(Subscript[z, o]),z[gtmin-t]]/.gminsol,(Subscript[z, f])},{t,-tplot/10,tplot},PlotRange->All,AxesLabel->{"",z},AxesOrigin->{-tplot/10,Automatic},DisplayFunction->Identity];
If[$VersionNumber==6.0,
Show[GraphicsGrid[{{mPlot,ExPlot,FPlot},{inPlot,knPlot,lnPlot},{itPlot,ktPlot,YPlot},{iPlot,χPlot,zPlot}}]],
Show[GraphicsArray[{{mPlot,ExPlot,FPlot},{inPlot,knPlot,lnPlot},{itPlot,ktPlot,YPlot},{iPlot,χPlot,zPlot}}]]]

Labels={"m","Ex","F","(Subscript[I, n])","(Subscript[k, n])","(Subscript[l, n])","(Subscript[I, t])","(Subscript[k, t])","Y","i","χ","z"};InitialSteadyState={(Subscript[m, o]),(Subscript[Ex, o]),(Subscript[F, o]),δ (Subscript[kn, o]),(Subscript[kn, o]),(Subscript[ln, o]),δ (Subscript[kt, o]),(Subscript[kt, o]),(Subscript[Y, o]),ρ+(Subscript[π, o]),(Subscript[π, o]),(Subscript[z, o])}//N;FinalSteadyState={(Subscript[m, f]),(Subscript[Ex, f]),(Subscript[F, f]),δ (Subscript[kn, f]),(Subscript[kn, f]),(Subscript[ln, f]),δ (Subscript[kt, f]),(Subscript[kt, f]),(Subscript[Y, f]),ρ+(Subscript[π, f]),(Subscript[π, f]),(Subscript[z, f])}//N;
(Dat=Table[{tt,m[t],Ex[Pn[t],λ1[t]],F[t],in[Pn[t],λ1[t],λ2[t],kn[t]],kn[t],ln[t],it[Pn[t],λ1[t],λ3[t],kt[t]],kt[t],Pn[t]F[kn[t],ln[t]]+z[t]G[kt[t],lbar-ln[t]],i[Pn[t],λ1[t],m[t],F[t]],χ[Pn[t],λ1[t],m[t],F[t]],z[t]}/.{t-> gtmin-t}/.tt->t/.gminsol,{t,0,31,1/5}])//MatrixForm;
(Dat=Join[{Join[{"Initial Steady State"},InitialSteadyState],Join[{"Final Steady State"},FinalSteadyState],Join[{Time},Labels]},Flatten[Dat,1]])//N//MatrixForm

The solve 4D Package  (* compare with 3D version and make changes *)

(* Interface *)
BeginPackage["solve`"];

EndPackage[];

Off[General::spell];Off[General::spell1];Off[Unset::norep];Off[FindMinimum::lstol];Off[FindMinimum::reged];Off[InterpolatingFunction::dmval];Off[Precision::precsm];

(* Usage messages for the exported functions and the context itself *)
{solve::usage = "solve[sys, {t}, vars, varso, varsf, varsendo, init] solves a system of first-order differential equations. \n
sys   = {xdot, ydot,...} lists the right hand side of the system of first-order differential equations (with state variables at end).\n
vars  = {x, y, ...} is the list of variables in the system in the same order as the differential equations in sys.\n
varso = {(Subscript[x, o]), (Subscript[y, o]), ...} is the list of original steady values of variables in the system in the same order as vars.\n
varsf = {(Subscript[x, f]), (Subscript[y, f]), ...} is the list of final steady values of variables in the system in the same order as vars.\n
varsendo = {{w,{x,y, ...}}, ...} is the list of the endogenous variables {w, ...} where w is function of {x,y, ...} where {x,y,...} are in the order in which they appear as arguments of w. If no endogenous variables enter {{,{}}}.\n
init   = {tol, precision, epsilon} sets the required accuracy to tol, starting precision to precision, and radius of initial ϵ-ring to epsilon." ;};

(* Error messages for the exported objects *)
{  solve::sysDESys = "You called `1` with first argument `2` (the system of differential equations) that has variables other than those defined in `3` and `4`.";

solve::sysLength = "You called `1` with lengths of `2`, `3`, `4` and `5` not equal.";

solve::vars = "You called `1` with third argument `2` that is not a list of variables.";

solve::origVals = "You called `1` with fourth argument `2` that is not a list of numbers.";

solve::finVals = "You called `1` with fifth argument `2` that is not a list of numbers.";

solve::endoVars = "You called `1` with sixth argument `2` that is not a list of lists of the form {w, {x,y, ...}} where w is an endogenous variable and {x,y, ...} are the variables in the system as defined in `3`.";

solve::initVals = "You called `1` with last argument `2` that is not a list of numbers.";

solve::sysError = "You called `1` with a system that does not have four negative eigenvalues when linearized."};

(* Implementation *)
Begin["solve`Private`"];

(* Definition of the exported functions *)
solve[asys_List,{at_},avars_List,avarso_List,avarsf_List,avarsendo___List,ainit_List]:=Block[{sys,t,vars,varso,varsf,varsendo,init,T,tol,precision,delta,Tmax,firstsim,firstrun,system,box,A,sysL,dir,localmin,pmin,epsilon,digits,ntrajini,nringsini,gmindis,gtmin,gminsol,initsol,ini,dis,tsol,tdis,doneflag},

{sys,{t},vars,varso,varsf,varsendo,init}={asys,{at},avars,avarso,avarsf,avarsendo,ainit};{T,tol,precision,delta}=init;Tmax=65;firstsim = 01;firstrun = 01;

(* Check input *)
checkInput[];

(* Make the system from input data and create box over which to solve the system *)
system=makeSys[];
box:=Apply[Or,Join[MapThread[First[Map[(#1),{t}]]<((Min[#2,#3])/10)&,{vars,varso,varsf}],MapThread[First[Map[(#1),{t}]]>2 Max[#2,#3]&,{vars,varso,varsf}]]];

(* Solve the linearized system *)
A=linearize[];
sysL=Thread[Map[ToExpression[StringJoin[ToString[#],"'","[",ToString[t],"]"]]&,vars]==-A.(Map[ToExpression[StringJoin[ToString[#],"[",ToString[t],"]"]]&,vars]-varsf),List];
dir = Sign[Take[vars,-4]-Take[varsf,-4]/.Thread[vars->sysLsol[t],List]/.t->1000//N];

(* SOLVE FOR THE GLOBAL NONLINEAR SADDLEPATH *);
While[True,
localmin=0;

Label[begin];
$MinPrecision=digits=precision;ntrajini=15;nringsini=7;gmindis=tol;

(* Finding the center and the radius of the initial ring for first run of search. *)
If[firstsim==1 && firstrun ==0,{pmin,epsilon}=(<<solution1)[[{3,4}]]];If[firstsim==0 && firstrun ==0, {pmin,epsilon}=(<<solution2)[[{3,4}]]]; 

(* Find the intersection of the linear solution with relevant face of the δ-diamond/simplex to use as the center. *)
If[firstrun ==1,
	epsilon=delta/20;
	
	initsol=FindRoot[dir.Take[sysLsol[t]-varsf,-4]==delta,{t,10},WorkingPrecision->digits];
	pmin=Take[sysLsol[t],-4]/.initsol;Print[{dir,initsol,pmin}];
	
	If[initsol[[1,2]] < T+10  && firstsim ==1,
		delta=delta/2;Print["The δ-diamond/simplex is too large. Decreasing 'delta'."];Goto[begin]];
	If[initsol[[1,2]] > T+20  && firstsim==1,
		delta = 3delta;Print["The δ-diamond/simplex is unnecessarily tight. Increasing 'delta' by a factor of 3."];Goto[begin]]; ];

(* Conduct an initial check for sufficiency of precision in calculations using the linearized system. *)
If[firstsim==1 && firstrun ==1,
	ini:=Thread[Table[vars[[i]][0],{i,1,Length[vars]}]==Apply[Sequence,{stabcond[pmin],pmin},2],List];firstit=1;box1=True;

	{dis,tdis,tsol}=compdisini[sysL,ini];
	
	Print["LINEARIZED SYSTEM: ","    Minimum Distance from Initial Steady State = ",dis//N,
		"    Time of Closest Approach = ",tdis//N,"    Solution Time = ",tsol//N];
	If[dis <1000 tol,
		Print["Preliminary check indicates you have enough precision for calculations."],precision=precision+4;
		Print["You may not have enough precision for calculations. Increasing 'precision' by 4."];Goto[begin]];];

(* Calculate the minimum distance for the trajectory of the nonlinear system that starts from initial point on the δ-diamond/simplex. *)
ini:=Thread[Table[vars[[i]][0],{i,1,Length[vars]}]==Apply[Sequence,{stabcond[pmin],pmin},2],List];firstit=1;box1=box;

{dis,tdis,tsol}=compdisini[system,ini];gmindis=dis;

Print["NONLINEAR SYSTEM: ","    Minimum Distance from Initial Steady State = ",dis//N,"    Time of Closest Approach = ",tdis//N,"    Solution Time = ",tsol//N];

(* SOLVE - Solve for the global nonlinear saddle path*) ;
If[localmin==1,pmin=pmin+delta/100{1/3,1/3,-2/3}];

{gtmin,gminsol,pmin,epsilon}=compdis3[pmin,epsilon];

If[gtmin==(10^100),
precision=precision+4;
If[precision<80,Print["Increased 'precision' by 4 digits to ",precision," digits."],
Print["Either the system is highly ill-conditioned or you might have encountered a local minima. Perturing the center of the initial ring and resetting 'precision' to its initial value"];precision=16;localmin=1];Goto[begin]];

(* ERROR CHECK - Check for the accuracy of the solution over the time period of interest. *) 
{firstsim,firstrun, delta, precision, doneflag}=errorCheck[];
If[doneflag==1,Goto[end]];

];
 
(* DURATION CHECK - Check the duration for which the solution is accurate within tolerance 'tol'.*);
Label[end];
durationCheck[];

{gtmin,gminsol}];

checkInput[___]:=Block[{},

If[Not[Length[sys]==Length[vars]==Length[varso]==Length[varsf]],
	Message[solve::sysLength,"'solve'","'sys'","'vars'","'varso'","'varsf'"];Abort[]];

If[Not[Length[vars]==Length[Cases[vars,_Symbol]]],Message[solve::vars,"'solve'","'vars'"];Abort[]];

If[Not[Length[varso]==Length[Cases[varso,_?NumericQ]]],Message[solve::origVals,"'solve'","'varso'"];Abort[]];

If[Not[Length[varsf]==Length[Cases[varsf,_?NumericQ]]],Message[solve::finVals,"'solve'","'varsf'"];Abort[]];

If[Not[Length[varsendo]==Length[Cases[varsendo,{_Symbol,_List}]]]||Not[Length[Flatten[varsendo]]==Length[Cases[Flatten[varsendo],_Symbol]]],
	Message[solve::endoVars,"'solve'","'varsendo'","'vars'"];Abort[];];

If[Not[Length[init]==Length[Cases[init,_?NumericQ]]],Message[solve::initVals,"'solve'","'init'"];Abort[]];

];

makeSys[___]:=Block[{system},

system=-sys/. Thread[MapThread[(#1)&,Transpose[varsendo]]->
					MapThread[ToExpression[StringJoin[ToString[ #1],"[",ToString[#2],"]"]] /.List->Sequence&, Transpose[varsendo]],
						List];

If[Length[Cases[system/.Thread[vars->varsf,List],_?NumberQ]]≠Length[vars],
	Message[solve::sysDESys,"'solve'","'sys'","'vars'","'varsendo'"];Abort[];];

system=Thread[Map[ToExpression[StringJoin[ToString[#],"'","[",ToString[t],"]"]]&,vars]==  
			(system=system/.Thread[vars->Map[ToExpression[StringJoin[ToString[#],"[",ToString[t],"]"]]&,vars],List]),
				List];

system];

linearize[___]:=Block[{ sysL,A,vals,vecs, negvalsPos, negvecsT,dev},

dev= (10^(-50));

sysL= sys/.Thread[MapThread[(#1)&,Transpose[varsendo]]->
			MapThread[ToExpression[StringJoin[ToString[ #1],"[",ToString[#2],"]"]] /.List->Sequence&, Transpose[varsendo]],
				List];
A=Transpose[Table[(sysL/.Thread[vars->varsf+Table[If[j==i,dev,0],{j,1,Length[vars]}],List])-(sysL/.Thread[vars->varsf+Table[If[j==i,-dev,0],{j,1,Length[vars]}],List]),{i,1,Length[vars]}]/(2dev)];
{vals,vecs}=Eigensystem[A];

negvalsPos=Flatten[Position[Map[Chop[Re[#]]<0&,vals],True]];
If[Length[negvalsPos]≠ 4,
	Message[solve::sysError,"'solve'"];
	Print["The eigenvalue of the linearized system are: "];
	Print[vals//N];Abort[]];
negvecsT = Transpose[Part[vecs,negvalsPos]];

sysLsol[x_]:=Evaluate[varsf-negvecsT.DiagonalMatrix[Map[Exp,Select[vals,Chop[Re[#]]<0&]x]].
						Inverse[Part[negvecsT,Table[i,{i,Length[vars]-3,Length[vars]}]]].Take[varsf-varso,-4]];

stabcond[{x_,y_,z_,w_}]:=Evaluate[First[Simplify[Drop[vars,-4]/.
						Solve[(Drop[vars-varsf,-4])==Part[negvecsT,Table[i,{i,1,Length[vars]-4}]].
						Inverse[Part[negvecsT,Table[i,{i,Length[vars]-3,Length[vars]}]]].({x,y,z,w}-Take[varsf,-4]),
							Drop[vars,-4] ,WorkingPrecision->Infinity]]]];
A];

compdisini[system_,ini_]:=Block[{sol,tsol,dis,tdis},

sol=NDSolve[{system,ini},vars,{t,0,Tmax},
		Method->{"EventLocator", "Event"->box1,"EventAction":> Throw[Null, "StopIntegration"],"EventLocationMethod"->"StepBegin"},
		MaxSteps->100000,WorkingPrecision->digits];

tsol=sol[[1,1,2,1,1,2]];{dis,tdis}=distance[sol,tsol];

{dis,tdis,tsol}];

distance[sol_,tsol_]:=If[firstit==1,
Block[{dis,dislast,tdis,steps,tbeg,tend,intsz,traj},

dis= (10^7);dislast=(10^8);steps=4 Floor[tsol] ;tbeg=0;tend=tsol;intsz=tend-tbeg;

While[ (dislast-dis)>gmindis (10^(-5)),
dislast=dis;
traj=Table[Norm[Flatten[{vars[[-4]][t]-varso[[-4]],vars[[-3]][t]-varso[[-3]],vars[[-2]][t]-varso[[-2]],vars[[-1]][t]-varso[[-1]]}/.sol],2],{t,tbeg,tend,(intsz/steps)}];
dis=Min[traj];tdis=tbeg+(Position[traj,dis][[1,1]]-1)(intsz/steps);
tbeg=Max[tdis-(intsz/steps),0];tend=Min[tdis+(intsz/steps),tsol];intsz=tend-tbeg;steps=10];
{dis,tdis}],

Block[{dis,tdis,tbeg,tend,min},

tbeg=0;tend=tsol;

min=FindMinimum[First[Norm[Flatten[{vars[[-4]][t]-varso[[-4]],vars[[-3]][t]-varso[[-3]],vars[[-2]][t]-varso[[-2]],vars[[-1]][t]-varso[[-1]]}]]/.sol],
				{t,9tsol/10,tbeg,tsol},WorkingPrecision->digits,AccuracyGoal->5-Log[10,gmindis]];

{dis,tdis}={min[[1]],min[[2,1,2]]}]];

(* Set up for solving for the global nonlinear saddle path *)
compdis3[p_,ϵ_]:=Block[{gmin,dis,tdis,sol,iflag,bit},

{pmin,epsilon}={p,ϵ};bit=0;

{pmin,gmindis,gtmin,epsilon}=compdis2[pmin,epsilon,-π,-(π/2)];

iflag=0;bit=1;

While[True, (*Label[begin];*)
Print[" "];Print["BEGINNING SEARCH ON BALL ",bit," WITH CENTER ",pmin];

{{dis,tdis,sol},gmin}=compdis1[pmin,epsilon,-π,-(π/2)];If[tdis==(10^100),gtmin=tdis;Break[];Goto[end]];

If[dis≤gmindis,gmindis=dis;gtmin=tdis;gminsol=sol;pmin=gmin;epsilon=(epsilon/1),epsilon=(epsilon/2)];
If[gtmin>Tss+10 && iflag==0,iflag=1;{gminsol,gtmin,pmin,5epsilon}>>solution1];
If[gtmin>Max[Tss+10 ,(<<solution1)[[2]]+5]&& iflag==1 ,
	{gminsol1,gtmin1,pmin1}=(<<solution1)[[{1,2,3}]];
	If[Max[sysdiff]>tol,
		Print[sysdiff];{gminsol,gtmin,pmin,5epsilon}>>solution1,
		{gminsol,gtmin,pmin,5epsilon}>>solution1;Print[sysdiff];Goto[end]];];

Print["AT THE END OF SEARCH ON BALL = ",bit,":    CURRENT MINIMUM DISTANCE = ",gmindis//N,"     TIME OF CLOSEST APPROACH = ",gtmin];bit++;

If[gmindis<tol/3,Print["Target Tolerance = ",tol/3,"   Tolerance Achieved = ",gmindis];Break[]; Goto[end]; ];(*Goto[begin];*)
];

Label[end];
{gtmin,gminsol,pmin,epsilon}];

compdis2[p_,ϵ_,thetabeg_,phibeg_]:=Block[{ibmin,ibtmin,ibminsol,ibpmin,ibit},

epsilon=ϵ;pini=p;trajtrack=0;ibit=1;

While[True,

{{ibmin,ibtmin,ibminsol},ibpmin}=compdis1[pini,epsilon,thetabeg,phibeg];
Print["INITIAL BALL = ",ibit,":    MINIMUM DISTANCE = ",ibmin//N,"     TIME OF CLOSEST APPROACH = ",ibtmin];ibit++;

(*Print[{{rmin,rtmin,rminsol},rpmin}];*)
If[ibmin≤gmindis,
	gmindis=ibmin;gtmin=ibtmin;pmin=ibpmin;
	If[firstrun==1,epsilon=epsilon/2,epsilon=2epsilon],
	If[firstrun==1,epsilon=(*4*) 1epsilon,epsilon=epsilon/4(*2*)]; Break[] ];
];
Print["AFTER INITIAL SEARCH:    MINIMUM DISTANCE = ",gmindis//N,"     TIME OF CLOSEST APPROACH = ",gtmin];

{pmin,gmindis,gtmin,epsilon}];

(* Set up for solving for the global nonlinear saddle path *)
compdis1[p_,epsilon_,thetabeg_,phibeg_]:=Block[{bpmin,bpmin1,bpmin2,bmindis,jb1,jb2,sol,dis,tdis,
btmin,jbmin,bminsol,btmin2,jbmin2,bminsol2,bmindis2,theta,jbdis,jbtdis,jbdissol,jbp,btab,bminpos,initring},

{nrings,bmindis,bmindis2}={nringsini,(10^7),(10^7)};trajtrack=0;initring=1;

(*Label[begin];*)
For[jb=1,jb≤ nrings ,jb++,
	theta= thetabeg+((2π jb)/nrings);
	
	{{dis,tdis,sol},bpmin}=compdis[p,epsilon,theta,phibeg];
	
	If[tdis==(10^100),gtmin=tdis;Goto[end]];
	jbdis[jb]=dis;jbtdis[jb]=tdis;jbdissol[jb]=sol;jbp[jb]=bpmin;
];


btab=Table[jbdis[Mod[ib-1,nrings]+1]-jbdis[Mod[ib,nrings]+1],{ib,1,nrings}];bminpos=Flatten[Position[Table[btab[[Mod[ib-2,nrings]+1]]>0 &&btab[[Mod[ib-1,nrings]+1]]<0,{ib,1,nrings}],True]];Print[" "];

For[ib=1,ib≤Length[bminpos],ib++,
	Print["Initial Ring = ",ib,"   {bmindis, btdis, bjmin, nrings} = ",
	{jbdis[bminpos[[ib]]],jbtdis[bminpos[[ib]]]//N,bminpos[[ib]],nrings}]
];Print[" "];initring=0;

For[ib=1,ib≤Length[bminpos],ib++,

{bmindis,btmin,jbmin,nrings,bminsol}={jbdis[bminpos[[ib]]],jbtdis[bminpos[[ib]]],bminpos[[ib]],nringsini,jbdissol[bminpos[[ib]]]};

{jb1,jb2}=2{jbmin-1,jbmin+1};{nrings,jbmin}=2{nrings,jbmin};bpmin1=jbp[bminpos[[ib]]];

Label[begin];
For[jb=jb1+1,jb≤ jb2  ,jb+=2,

theta= thetabeg+2π jb/nrings;{{dis,tdis,sol},bpmin}=compdis[p,epsilon,theta,phibeg];

If[tdis==(10^100),gtmin=tdis;Goto[end]];
If[ Abs[bmindis-dis]≤Max[(10^(-6)),(gmindis/100)],Goto[end] ];
If[dis≤bmindis,bmindis=dis;btmin=tdis;jbmin=jb;bminsol=sol;bpmin1=bpmin;Break[];];
];

{jb1,jb2}=2{jbmin-1,jbmin+1};{nrings,jbmin}=2{nrings,jbmin}; 
Goto[begin];

Label[end];
If[bmindis<bmindis2,
	bmindis2=bmindis;btmin2=btmin;jbmin2=jbmin;bminsol2=bminsol;bpmin2=bpmin1];
If[bmindis2<gmindis,Break[]];(* UNDER TEST: to test if on finding improvement on a ball while searchin one initial ring you can skip searching on the other initial rings. *)
];

If[bit≠0,Print["Ball = ",bit,":    Minimum Distance = ",bmindis2//N,"    Time of Closest Approach = ",btmin2//N,"    No. of Subdivisions of ϵ-ball > ",jbmin2];];

{{bmindis2,btmin2,bminsol2},bpmin2}];

(* Set up for solving for the global nonlinear saddle path *)
compdis[p_,ϵ_,theta_,phibeg_]:=Block[{pini,pmin1,pmin2,mindis,j1,j2,firstit,sol,tsol,dis,tdis,tmin,jmin,minsol,tmin2,jmin2,minsol2,mindis2,jdis,jtdis,jdissol,jp,solveit,ntraj},

mindis=mindis2=(10^7);j1=1;j2=ntraj=ntrajini;firstit=1;

solveit:=(phi=phibeg+((1π j)/ntraj);pini=p+ϵ(Sin[phi]dir vec[[1]]+ Cos[theta]Cos[phi]dir vec[[2]]+Sin[theta]Cos[phi]dir vec[[3]]);
ini=Thread[Table[vars[[i]][0],{i,1,Length[vars]}]==Apply[Sequence,{stabcond[pini],pini},2],List];
sol=NDSolve[SetPrecision[{system,ini},precision],vars,{t,0,Tmax},Method->{"EventLocator", "Event"->SetPrecision[box,precision],"EventAction":> Throw[Null, "StopIntegration"],
"EventLocationMethod"->"StepBegin"},MaxSteps->10000,WorkingPrecision->digits];{dis,tdis}=distance[sol,sol[[1,1,2,1,1,2]]];);

For[j=1,j≤ ntraj  ,j++,
	If[trajtrack==0||(trajtrack==1 && (j==minpos[[1]]||j==minpos[[2]])),
		solveit;jdis[j]=dis;jtdis[j]=tdis;dissol[j]=sol;jp[j]=pini ,jdis[j]=(10^7)];
];

tab=Table[jdis[ii]-jdis[ii+1],{ii,1,ntraj-1}];firstit=0;
minpos=Flatten[Position[Table[If[ii==1,tab[[ii]]<0,If[ii==ntraj,tab[[ii-1]]>0,tab[[ii-1]]>0 &&tab[[ii]]<0]],{ii,1,ntraj}],True]];

If[rit≠0 &&Length[minpos]==2,trajtrack=0];
(*If[(*rit≠ 0*)True,
Print["Local minima from initial search over the ring."];
For[ii=1,ii≤Length[minpos],ii++,
Print[ii,":" ,"   {mindis, tdis, jmin, ntraj} = ",{jdis[minpos[[ii]]],jtdis[minpos[[ii]]]//N,minpos[[ii]],ntraj}]
];];*)

For[ii=1,ii≤Length[minpos],ii++,

{mindis,tmin,jmin,ntraj,minsol}={jdis[minpos[[ii]]],jtdis[minpos[[ii]]],minpos[[ii]],ntrajini,dissol[minpos[[ii]]]};
pmin1=jp[minpos[[ii]]];{j1,j2}=2{jmin-1,jmin+1};{ntraj,jmin}=2{ntraj,jmin};

While[True,(*Label[begin];*)
For[j=j1+1,j≤ j2  ,j+=2,
	If[ ntraj>(10^(precision+15)),Print["You may be losing precision in calculations."] ;tmin2=(10^100);Goto[end2];];
	solveit;
	If[ Abs[mindis-dis]≤Max[(10^(-8)),(gmindis/10000)],Goto[end1] ];
	If[dis≤mindis,
		If[phi<phibeg||phi>phibeg+π,mindis=(10^7);Goto[end1]];
		mindis=dis;tmin=tdis;jmin=j;minsol=sol;pmin1=pini;Break[]];];

{j1,j2}=2{jmin-1,jmin+1};{ntraj,jmin}=2{ntraj,jmin};
 (*Goto[begin];*)
];

Label[end1];
If[mindis<mindis2,
	mindis2=mindis;tmin2=tmin;jmin2=jmin;minsol2=minsol;pmin2=pmin1];
If[initring==0 && mindis2<gmindis,Break[]];(* UNDER TEST : to check if it is okay to break off from ring search on finding improvement 
												if it is not the initial ring of the ball. initring is set in compdis1. *);
];

(*Print[pmin2];*)
If[(*rit≠ 0*)True ,
	Print["Ring = ",jb,"    Minimum Distance = ",mindis2//N,"    Time of Closest Approach = ",tmin2//N,"    No. of Subdivisions of ϵ-ring > ",jmin2];];

Label[end2];
{{mindis2,tmin2,minsol2},pmin2}];

vec={{-0.4872208210285517339656346783519891269808160803797712415865039589702486887293003407807400845236750724596019810623729`100.,-0.4999394618730551674197088470317218156155844750866232835905547349469924022821306227294551758926722320533613913660236`100.,0.3808106089032299867888231569346050257196958282990104583876035386229002093695948320656574474515966807631674044171052`100.,0.6063496739983768035742179125826457403700801217957550794927550587230441480043964164358161978761160337023595431787359`100.},{-0.7151123073705422758984923213248647968590500088695194130261061755424791564857436296721273488694944658068633868667997`100.,0.6962749513563497998843390466911133690677024742322083713736336137378054594826899203735843489048761502721949862665217`100.,0.0520721011201536841367114678382023907307452980744588325879757555920404974484140610227323366981588865981859982124163`100.,-0.0332347451059610277113167024214531761161327797082506865783655368590096082845208148350167122145093622364334072592356`100.},{0.0350795461563850278118692916514201197719367779270019041801732213263542380175013983498226088466672483779825055343508`100.,0.1235396362747044612207540037487089999635319873888117543866955503276619017682651950901608951264033722373964568194175`100.,-0.7760619122124071480872482322238330899084646356724542942162567299548466901260571052297105295241890777688061890323372`100.,0.617442729781317090065324850508095565576544767827042075753799963372934790448411972370028748221866183160315549411722`100.}};

errorCheck[___]:=Block[{diffprec,diffdel, satisprec, notsatisprec, notsatisdel, increaseprec,decreasedel,
step2, step3, congratulations, done,gtmin1,gminsol1},

doneflag=0;

(* Messages *)
diffprec:=Print["Difference in values of variables due to latest increase in precision = ", sysDiff[0,T]//N];
diffdel:=Print["Difference in values of variables due to latest decrease in delta = ",sysDiff[0,T]//N];

satisprec= "THE PRECISION IS SATISFACTORY. The variables changed by less than 'tol' due to latest increase in 'precision' for current δ-diamond\\simplex.";
notsatisprec="The precision is NOT satisfactory. The variables changed by greater than 'tol' due to latest increase in 'precision' for current δ-diamond\\simplex.";
increaseprec:=Print["Increased 'precision' by 4 digits to ",precision," digits."];

notsatisdel="The δ-diamonds are NOT tight enough. The variables changed by greater than 'tol' due to latest decrease in 'delta'.";
decreasedel:=Print["Decreased 'delta' by a half to ",N[delta,precision],"."];

step2="Starting Step 2: Checking for the precision of the candidate solution.";step3="Starting Step 3: Checking for the accuracy of the candidate solution.";

congratulations= "CONGRATULATIONS! The variables changed by less than 'tol' due to latest decrease in 'delta'.";
done = "You are done! Hope to see you in another model. Continue with the next cell.";

(* Check candidate solution *)
Which[
firstsim == 1 && firstrun == 1, 
	{gminsol,gtmin,pmin,5epsilon}>>solution1;Print[step2];precision=precision+4;increaseprec;firstrun=0,

firstsim == 1 && firstrun == 0,
	{gminsol1,gtmin1}=(<<solution1)[[{1,2}]];{gminsol,gtmin,pmin,5epsilon}>>solution1;diffprec;
	If[Max[sysDiff[0,T]]>tol,Print[notsatisprec];precision=precision+4;increaseprec,Print[satisprec];Print[step3];
		firstsim=0;firstrun=1;delta=delta/2;decreasedel],

firstsim == 0  && firstrun == 1,
	{gminsol,gtmin,pmin,5epsilon}>>solution2;{gminsol1,gtmin1}=(<<solution1)[[{1,2}]];diffprec;
	If[Max[sysDiff[0,T]]>tol,firstrun=0;precision=precision+4;increaseprec,diffdel;Print[congratulations];
		Print[done];doneflag=1],

firstsim == 0  && firstrun == 0,
	{gminsol1,gtmin1}=(<<solution2)[[{1,2}]];{gminsol,gtmin,pmin,5epsilon}>>solution2;diffprec;
	If[Max[sysDiff[0,T]]>tol,
		Print[notsatisprec];precision=precision+4;increaseprec,
		Print[satisprec];{gminsol1,gtmin1}=(<<solution1)[[{1,2}]];diffdel;
		If[Max[sysDiff[0,T]]>tol,
			{gminsol,gtmin,pmin,5epsilon}>>solution1;Print[notsatisdel];firstrun=1;delta=delta/2;decreasedel,
			Print[congratulations];Print[done];doneflag=1];];
];   

{firstsim,firstrun,delta,precision,doneflag} ];

sysDiff[t1_?NumberQ,t2_?NumberQ]:=Table[Max[Abs[Table[((vars[[i]][gtmin-t]/.gminsol)-(vars[[i]][gtmin1-t]/.gminsol1))[[1]],{t,t1,t2,1/10}]]],{i,1,Length[vars]}];

durationCheck[___]:=Block[{gminsol1,gtmin1,j},

{gminsol1,gtmin1}=(<<solution1)[[{1,2}]]; {gminsol,gtmin}=(<<solution2)[[{1,2}]];

For[j=T+1,j≤IntegerPart[gtmin1],j++,
	If[Max[sysDiff[j-1,j]]>tol,
		Print["The solution is accurate within tolerance 'tol' upto ",j-1," years. Continue with next cell."];Break[]];];

If[j==IntegerPart[gtmin1]+1,Print["The solution is accurate upto ",j-1," years. Continue with next cell."]];

];

End[];

(* Epilog *)
Protect[Evaluate[$Context<>"*"]];